clear
close all
clc

%% Parametros de ajuste de la simulacion
M = 32;  % cantidad de simbolos
k = log2(M); % bits por simbolo
SNRbdB = 2:14; % SNR por bit en dB
nsimbs = 4e5; % cantidad de simbolos que para cada valor de SNRbdB en el Montercarlo


%% Conversiones e inicializaciones varias
SNRbveces = 10.^(SNRbdB/10); % SNR por bit en veces
SNRdB = SNRbdB + 10 * log10(k);  % SNR en dB
SNR = 10.^(SNRdB/10);    % SNR en veces
Eavg = 1;

Pes = zeros(size(SNRbdB));  % probabilidad de error por simbolo simulada
Peb = zeros(size(SNRbdB));  % probabilidad de error de bit simulada


%% P(error) teorica:
Pet = @(SNR) 4 * (1-1/sqrt(M)) * qfunc(sqrt(3/(M-1)/2*SNR)) .* (1- (1-1/sqrt(M)) * qfunc(sqrt(3/(M-1)/2*SNR)));

%% Simulacion Montecarlo

% por su energia promedio Eavg (definida mas arriba)
modPAM = comm.RectangularQAMModulator(M,'SymbolMapping','Gray',... % codifica en Gray %
    'NormalizationMethod','Average Power',...
    'AveragePower',Eavg); % Energia promedio deseada = 1

% Genero un objeto "demodulador PAM" con las mismas especificaciones que el
% modulador
demodPAM = comm.RectangularQAMDemodulator(M,'SymbolMapping','Gray',... % codifica en Gray 
    'NormalizationMethod','Average Power',...
    'AveragePower',Eavg); % Energia promedio deseada


for n0 = 1:length(SNRbdB)
    
    %% Primero trabajo con simbolos y despues analizo por bit
    
    sigma = sqrt(Eavg/SNR(n0)); % desvio de ruido
    simbtx = randi(M,nsimbs,1)-1;   % simbolos a transmitir (en 1:M);
    
    % la funcion "step" llama al objeto modulador o demodulador
    simbmod = step(modPAM, simbtx); % simbolos modulados para su transmision "OJO, son COMPLEJOS"
      
    % Sumo el ruido (OJO sumo el ruido como si la constelacion tuviera parte compleja)
    simbmod = simbmod + sigma * (randn(nsimbs,1)+1i * randn(nsimbs,1));
    
    % Demodulo
    simbrx= step(demodPAM,simbmod); % simbolos recibidos
    
    % Probabilidad de error de simbolo simulada
    Pes(n0) = nnz(simbtx~=simbrx)/nsimbs;
    
    %% Ahora analizo los errores de bit (el Gray lo hizo el modulador)
    
    bitstx = de2bi(simbtx); % convierto los simbolos transmitidos a bits
    bitsrx = de2bi(simbrx); % convierto los simbolos recibidos a bits
    % comparo los bits y cuento cuantos difieren
    Peb(n0) = nnz(xor(bitstx,bitsrx))/k/nsimbs; 
    
end


%% Graficos
% Probabilidad de error de simbolo
figure(1)
semilogy(SNRbdB,Pet(SNR),'-dr','linewidth',1.5);
hold on
semilogy(SNRbdB,Pes,'-o','linewidth',1.5);
grid on
box on
xlabel('SNR por bit (dB)');
ylabel('P_e')
legend('Curva teorica','Montecarlo')

% Probabilidad de error por bit
figure(2)
semilogy(SNRbdB,Pet(SNR)/k,'-dr');
hold on
semilogy(SNRbdB,Peb,'-o','linewidth',1.5);
grid on
box on
xlabel('SNR por bit (dB)');
ylabel('P_b')
legend('Curva teorica','Montecarlo')

